FOXSI Optics Simulation

by Steven Christe



In [1]:

    
from foxsisim.module import Module
from foxsisim.detector import Detector
from foxsisim.source import Source
from foxsisim.plotting import plot



In [2]:

    
import matplotlib.pyplot as plt
%matplotlib inline



In [4]:

    
import foxsisim

Define a point x-ray source located at a distance z in front of the optics



In [5]:

    
source_distance = 1e4
source = Source(type='point', center = [0, 0, -source_distance], color = [1,0,0])

Now define an energy spectrum for the source (this definition is a bit strange because of a "bug" in piecewise)



In [6]:

    
max_energy = 20.0
def spectrum(z):
        if (type(z) is not type([1])) and (type(z) is not type(np.array(1))):
            x = np.array([z])
        else:
            x = np.array(z)
        return np.piecewise(x, [x < 0, (x < max_energy) & (x > 0), (x >= max_energy)], [0, 1 / np.float(max_energy), 0])



In [7]:

    
source.loadSpectrum(spectrum)

Let's plot the spectrum



In [8]:

    
plot(source)

Now define a FOXSI telescope module



In [11]:

    
radii = [5.15100,4.90000,4.65900,4.42900,4.21000,4.00000,3.79900] # 7 shell radii
seglen = 30
base = [0,0,0]
focal_length = 200
module = Module(radii=radii, seglen=seglen, base=base, angles=None, focal=focal_length)









    



---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-11-c6f6c44aa3aa> in <module>()
      3 base = [0,0,0]
      4 focal_length = 200
----> 5 module = Module(radii=radii, seglen=seglen, base=base, angles=None, focal=focal_length)

/Users/schriste/Dropbox/Developer/python/foxsi-optics-sim/src/foxsisim/module.py in __init__(self, base, seglen, focal, radii, angles)
     44         for i, r in enumerate(radii):
     45             self.shells.append(Shell(base=base, focal=focal, seglen=seglen, ang=angles[i],
---> 46                                      r=r))
     47 
     48         # inner core (blocks rays going through center of module)

TypeError: __init__() got an unexpected keyword argument 'focal'

A plot of the cross-section of the telescope module



In [38]:

    
# generate cross section
fig1 = plt.figure(figsize=(9,3))
axes1 = fig1.gca()
module.plot2D(axes1,'b')
plt.vlines(seglen,-10,10, label='z=0', linestyles='dashed')
plt.vlines(focal_length+seglen,-10,10, label='focus', linestyles='dashed')
plt.legend()









    Out[38]:





<matplotlib.legend.Legend at 0x10f3efc10>

Now let's define the detector that will detect the rays cast through the telescope module



In [39]:

    
# create a detector located slightly in front of focal point
detector = Detector(center=[0,0,230]) # focal point is at 230 by default

The following runs the simulation by generating rays at the source and casting them through the optic. Finally the rays are caught by the detector.



In [56]:

    
nrays = 1000



In [57]:

    
#  generate nrays from the source
rays = source.generateRays(module.targetFront, nrays)



In [42]:

    
#  generate nrays from the source
rays = source.generateRays(module.targetFront, nrays)

Let's compare the spectrum of the generated rays to the source spectrum



In [48]:

    
energies = [ray.energy for ray in rays]



In [49]:

    
plt.figure()
plt.hist(energies)
plt.title('Input Ray Spectrum N = ' + str(len(energies)))
plt.show()



In [50]:

    
# pass rays through module
module.passRays(rays, robust=True)



In [17]:

    
# pass rays through module
module.passRays(rays, robust=True)



In [51]:

    
# catch rays at detector
detector.catchRays(rays)



In [19]:

    
# catch rays at detector
detector.catchRays(rays)



In [54]:

    
fig = plt.figure(figsize=(5,5))
plot(detector, energy_range=[0,20])









    












    





<matplotlib.figure.Figure at 0x10f3ea910>

A plot of the spectrum as seen by the detector



In [55]:

    
fig = plt.figure(figsize=(5,5))
axes = fig.gca()
detector.plotSpectrum(axes, range=[0,20], dE=1)



In [25]:

    
def simulate_foxsi(source_distance=1e6, nrays=500):

    source = Source(type='point', center = [0, 0, -source_distance], color = [1,0,0])
    source.loadSpectrum(spectrum)

    
    #  generate nrays from the source
    rays = source.generateRays(module.targetFront, nrays)

    # pass rays through module
    module.passRays(rays, robust=True)

    # catch rays at detector
    detector.catchRays(rays)

    return detector



In [5]:

    
import foxsisim.reflectivity as ref
r = ref.Reflectivity()
plot(r)
plt.show()



In [13]:

    
energy_range = r.energy_range()
angle_range = r.angle_range()
energies, angles = np.mgrid[energy_range[0]:energy_range[1]:100j,
                                    angle_range[0]:angle_range[1]:200j]
z = r.value(energies, angles)
extent = (energy_range[0], energy_range[1],
        angle_range[0] + r._points[1, 1], angle_range[1])



In [25]:

    
plt.figure()
plt.yscale('log')
plt.xscale('log')
levels = np.arange(0, 1, 0.1)
cs = plt.contour(z, levels, linewidths=2, extent=[0.1,20,1,3])
plt.clabel(cs, inline=1)
plt.ylabel('Energy [keV')
plt.xlabel('Angle [deg]')

energy = np.arange(0.1,20, 0.1)
angles = np.arange(0.1, 2, 0.1)

plt.figure()
angles = np.arange(1,200,1)
for this_angle in angles:
    plt.plot(z[:,this_angle])
plt.show()



In [20]:

    
plt.figure()
energy = np.arange(0.1,20, 0.1)
angles = np.arange(0.1, 2, 0.1)
for this_angle in angles:
    plt.plot(energy, r.value(energy,this_angle), label=str(this_angle))
plt.legend()
plt.show()



In [26]:

    
r.energy_range()









    Out[26]:





array([  1.,  30.])



In [27]:

    
r.angle_range()









    Out[27]:





array([  0.,  10.])



In [28]:

    
r._values









    Out[28]:





array([  1.00000000e+00,   9.91481000e-01,   9.83034000e-01, ...,
         1.08969000e-09,   1.08049600e-09,   1.07231100e-09])



In [29]:

    
r.material









    Out[29]:





'Ni'



In [ ]: